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Spurious instabilities in multiangle simulations of collective flavor conversion 
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The dense neutrino flux streaming from a core-collapse supernova can undergo self-induced flavor 
conversion caused by neutrino-neutrino refraction. Numerical studies of these nonlinear effects are 
challenging because representing the neutrino radiation field by discrete energy and angle bins can 
easily lead to unphysical solutions. In particular, if the number of angle bins Na is too small, flavor 
conversion begins too deep and produces completely spurious results. At the same time, Na = 1 
(single-angle approximation) can be a good proxy for the A'^a oo limit. Based on a linearized 
stability analysis, we explain some of the puzzling effects of discrete angle distributions. 

PACS numbers: 14.60.Pq, 97.60.Bw 



I. INTRODUCTION 

The almost freely streaming flux of neutrinos emitted 
from a collapsed supernova (SN) core starts at a radius 
of 10-30 km, depending on the explosion phase and the 
neutrino energies, and at a density of 10^^"^'^ g cm~'^. In 
such conditions, neutrino refraction in matter is so large 
that the propagation eigenstates are almost identical to 
the weak interaction eigenstates, even though all neutrino 
mixing angles are relatively large. Despite this simple 
boundary condition, the subsequent neutrino flavor evo- 
lution is surprisingly complicated due to the nonlinear 
impact of neutrino- neutrino refraction P, Q . 

One major problem in understanding this phenomenon 
is that numerical studies are challenging. If the integra- 
tion begins close to the neutrino-emitting region of the 
SN (the "neutrino sphere" ) , large oscillation frequencies 
caused by the matter effect require many time steps. 
Moreover, the neutrino radiation field must be repre- 
sented in some numerical form, in practice a finite num- 
ber of energy and angle bins. Thus far axial symmetry of 
the neutrino stream was always assumed, implying that 
there was only one angle variable. Even with this simpli- 
fication, reliable results require a large-scale computing 
effort [|. 

Far away from the SN core, for example for a detector 
on Earth, the angle distribution is irrelevant and good en- 
ergy resolution is sufficient. In particular, one can study 
the sharp spectral features that self-induced flavor con- 
version can produce [l|, 043 • Therefore, many studies 
used the single- angle approximation, i.e., all neutrinos 
are put in a single angle bin corresponding to emission 
at the neutrino sphere at one specific angle relative to 
the radial direction. The simplest assumption of radially 
emitted neutrinos is not possible because parallel-moving 
neutrinos do not provide a refractive effect on each other. 
We always represent the single-angle case by assuming 
neutrino emission at 45° at the neutrino sphere. 

To go beyond single-angle studies, upgrading to a small 
number Na of angular bins is not enough. Beginning 
with the earliest studies [l[, all multiangle simulations 
were haunted by the same peculiar effect. Although sim- 
ulations with a large number of angular bins (Na — >■ oo) 



and those with just a single angle {Na — 1) often provide 
similar results, simulations using a relatively small Na 
lead to results completely different from the two limiting 
cases. The required Na to avoid spurious solutions had 
completely eluded explanation. Depending on circum- 
stances, a few tens of modes might be enough, whereas 
in other cases one needs thousands or more bins. 

For a simple ("single-crossed" Q) neutrino spectrum, 
which is relevant during the SN accretion phase, self- 
induced flavor conversion begins at a critical onset radius. 
If only neutrino- neutrino refraction is relevant, then at 
larger neutrino densities, close to the SN core, the system 
is stable ("sleeping top" phase ^, |9(]). The onset radius 
turns out to be almost the same in the single-angle and 
many-angle limits. On the other hand, if Na > 1 is 
too small, flavor conversion begins at a smaller radius 
(where higher neutrino densities dominate) and tends to 
cause kinematical decoherence among angle modes [l0| . 
Similar spurious solutions appear also when the normal 
matter effect is important. 

This situation is to be contrasted with the role of en- 
ergy bins that can be chosen to the desired resolution 
without affecting the qualitative behavior. Notice that 
neutrino energies (or rather the vacuum oscillation fre- 
quencies) appear in the linear part of the Hamiltonian, 
responsible for the vacuum oscillations. The angle vari- 
ables, on the other hand, appear in the neutrino-neutrino 
interaction part which is the source of nonlinear effects. 
In the single-angle multienergy case, the flavor evolution 
is described by a small number of collective variables 
("A^-mode coherence"), independently of the number of 
energy bins [TTI - [l3j . The ultimate source for this simple 
behavior is the simplicity of the single-angle Hamilto- 
nian which contains as many constants of the motion as 
variables and thus is integrable [l^, ■ If detailed en- 
ergy resolution is not important, we can therefore study 
conceptual aspects of multiangle effects in the monochro- 
matic approximation where we use only one energy bin 
for neutrinos and one for antineutrinos. 

The observed flavor conversion effect in the case of few 
angle modes must begin with an instability, i.e., with a 
runaway flavor mode in the interacting neutrino system. 
Understanding self-induced flavor conversion based on a 
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linearized stability analysis was pioneered by Sawyer 15 1 
and further developed by our group and collaborators [16, 
[iTj . We here apply this technique to the 'Wq effect" and 
find that many puzzling numerical observations easily fall 
into place. The spectrum of runaway modes is indeed 
very different for 'Wa = few" from Na ~ 1 and Na — > 
oo. Additionally, the two neutrino mass hierarchies show 
striking differences. 



II. LINEARIZED STABILITY APPROACH 
A. Equations of motion 

Following earlier works [H, [13, [ll], we describe the 
two-flavor neutrino field by energy- and angle-dependent 
2x2 matrices $£;.„(r). Boldface characters denote matri- 
ces in flavor space. The diagonal ^e,u elements are the 
ordinary number fluxes ^ (flavor a) integrated over 
a sphere of radius r, with negative E for antineutrinos. 
Moreover, we use the "flavor isospin convention" where 
^E.u < antineutrinos {E < 0) and F^ ,^ > for 
neutrinos {E > 0). 

The off-diagonal elements of $E.„(r), which are ini- 
tially very small, represent phase information caused by 
flavor oscillations. The flavor evolution follows from the 
Schrodinger-like equation [l^ 



(1) 



(2) 



with the Hamiltonian matrix 

1 /]VT2 

ilE,u - — — + V2GfN, 
Vu \ 2E 

+ -; — ^ / dE' / du' ■ 

47rH Jo VuVu- 

The matrix of neutrino mass squares causes vacuum 
flavor oscillations and that of net charged lepton den- 
sities — diag(7ie— ng, n^— rt^, tIt-— ) adds the mat- 
ter effect. The third term provides neutrino- neutrino 
refraction. A neutrino radial velocity at radius r is 
w„ = (1 — uT(^ jr^Y^"^ , where R is the radius of the 
neutrino sphere, at which we label the neutrino angle 
modes by their emission angle dp.- Our angle variable u 
is then defined by Vu\r=R = cosOr = (1 — u)^/"^, equiv- 
alent to u = sin^ 6fi. The factor 1 — WtiWt,' comes from 
the current-current nature of the weak interaction and 
leads to multiangle effects. Moreover, Vu appears in the 
denominator because we follow the flavor evolution pro- 
jected on the radial direction. 

Next, we study the instability driven by the atmo- 
spheric Ato^ and the mixing angle ^13. In the relevant 
SN region, propagation eigenstates are almost identical 
with weak-interaction eigenstates unless self-induced fla- 
vor conversion occurs. This means that initially the off- 
diagonal elements of ^e.u are very small, justifying a lin- 
earized stability analysis, but otherwise we do not need 
a specific numerical value of the mixing angle. 



We switch to the variable w = Am^/2_E which is much 
better suited to study flavor oscillation than energy itself. 
Finally, the flux matrices are written in the form 



Tr*„ 



s: 



(3) 



where „ and „ are the flavor fluxes at the neu- 
trino sphere. The flux summed over all flavors, Tr$^_,j, 
is conserved in our free-streaming limit. The Ve survival 
probability is ^[1 + Stj.„(r)] in terms of the "swap fac- 
tor" — 1 < Suj^uir) < 1. The off-diagonal element S^.u, is 
complex and ^, -I- |>S'i^,mP — 1. 

We have formulated our equations so that, for Ato^ > 
0, they describe the inverted hierarchy (IH) of the two 
possible neutrino mass orderings. The normal hierar- 
chy (NH) can be implemented with the substitution 
Am^ — >■ — Am^ or, equivalently, w — > —lo. 



B. Stability condition 

The possible onset of self-induced flavor conversions 
is best described in terms of the complex numbers Suj,u 
which are very small as long as neutrinos are in the eigen- 
states of the weak interaction in the presence of matter. 
The small-amplitude limit means IS*!^,!!! ^ 1 and to lin- 



ear order 



1 . Assuming in addition a large distance 



from the source so that 1—Vu ^ 1, the evolution equation 
linearized in Suj.u and in u is |l^] 

i^r-S'w.u = + uX) S^,u 

— /i / du du! (u + u ) g^'u' Su^'.u' ■ (4) 



Here, .„ is the neutrino spectrum (a; < for antineu- 
trinos) which we normalize to the antineutrino flux, i.e. 
J_^duj Jg dug^^u — ~1- The "asymmetry" between 
neutrinos and antineutrinos is e = J duj du g^j.u- 

Refractive effects are encoded in the r-dependent pa- 
rameters 



A 



V2Gf [n,{r)-n,{r)]^, 



V2Gf[F,AR)-FM] 

47rr2 2r2 ' 



7r?l, (5) 



and often combined in A = A -I- e/i. The factor R^/2r^ 
means that only the multiangle impact of the neutrino- 
neutrino and matter effects are relevant for our stabil- 
ity analysis, not the densities themselves. We normalize 
the neutrino-neutrino interaction strength fi, and conse- 
quently the spectrum g^i^u, to the Vf,-Vx flux difference at 
a chosen radius R, the nominal neutrino sphere. Physical 
results do not depend on the choice of R. 

Writing solutions of the linear differential equation, 
Eq. (H]), in the form 



5. 



(6) 
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with complex frequency = 7 + Ik and eigenvector Q^^u 
leads to the eigenvalue equation [l^ 

{uj+uX-Vt) QuiM ^ [I \ dv! duo' {u+u') g^iu' Qu'.u- ■ (7) 



The right-hand side of this equation is a linear polyno- 
mial in u so that the eigenvector must have the form 



bu 



ui + uX — rt ^ 



(8) 



where a and b are complex numbers. The form of the 
eigenfunctions is that of a Mobius transformation in 
the complex plane. This means that for fixed u we 
have a circle parametrized by ui and for fixed uj a cir- 
cle parametrized by u. The physical range u G (0, 1) 
then maps to a circular arc in the complex plane. 

Following Ref. [l^ we note that, after inserting Eq. ^ 
into Eq. ([7]), both sides are linear polynomials in u. Self- 
consistency requires 



where 



/i - 1 h 



fi / duduj ■ 



u 9u 



UJ - 



(9) 



(10) 



In contrast with Ref. [l6| we have included a factor /i in 
the integral to make In dimensionless. 

Nontrivial solutions for a and b exist if the determinant 
of the matrix in Eq. Q vanishes, implying 



(/i - If = loh 



(11) 



Changing the neutrino mass hierarchy from inverted to 
normal is simply achieved by the sign change cj — > — cj in 
the denominator of the integrand of Eq. ((TOl) . 



isotropically into space without limb darkening. This as- 
sumption corresponds to a box spectrum in the u variable 
of the form 



Biu) = 



1 for < li < 1, 



(12) 



otherwise. 
Therefore, overall we study the neutrino spectrum 
9u,u = [-S{-uJo - w) + (1 + e) 6{ujo - uj)] B{u) . (13) 

B. Inverted-hierarchy solution 

To solve for the eigenvalues with the spectrum 
Eq. (1131) one can perform the integrals In analytically, 
leading to expressions involving logarithms and the arc- 
tan function. It is then straightforward to solve Eq. pT|) 
for the eigenvalues fl numerically. In the absence of mat- 
ter (A = 0) and for e = 1/2 we find the "Continuum (IH)" 
growth rate k = Im(r2) shown in Fig. [1] The interacting 
neutrino stream is stable for fi above a critical value fj,2 
and below another /xi. 

It is instructive to compare this continuum result with 
the single-angle approximation where we replace the 
box spectrum with a single mode placed at its center: 
B{u) — > (5(1/2 — u). This example corresponds to the 
original "flavor pendulum" Q. One can solve the eigen- 
value equation explicitly and finds 



yi-(2 + e)M+(f (14) 



The solution has an imaginary part for fj, between the 
values /zi,2 = (1 + e/2 ± yT+iy^ . For e = 1/2 we find 
/xi_2 = 20 ± 4V24. The maximum growth rate and the 
corresponding /x value is 



2y/T 



e _2(2 + e) 

— at /imax — 5 ■ (1") 



III. CONTINUUM VERSUS SINGLE ANGLE 
A. Neutrino spectrum 

We use the simplest nontrivial setup that allows us 
to study the role of discrete angle modes, i.e., we con- 
sider a neutrino fiux streaming from a sphere at radius 
R that emits monochromatic fiuxes of and Ve- The 
Vf, flux is taken to be 1 -|- e times the 9^ flux, represent- 
ing deleptonization. The chosen vacuum frequency for 
monochromatic neutrinos and antineutrinos ibwo deter- 
mines the frequency scale of the system. We simplify the 
calculations by choosing wq = 1 as the unit of measure 
for all other frequencies such as /i. A, k and 7. 

The angle distribution is taken to be black-body- 
like; i.e., the neutrino sphere is taken to emit neutrinos 



K 




FIG. 1: Growth rate k as a function of effective neutrino 
density jj, in the absence of matter (A = 0) and assuming 
e — 1/2. We show the single-angle and continuum cases for 
normal and inverted hierarchy. 
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For e = 1/2 we find Kmax = v24 and /Xmax = 20. In other 
words, while wq is the only frequency scale in our problem 
and the dimensionless parameter e is not especially small, 
the /i range where the system is unstable reaches to fi2 ^ 
cjQ and likewise, typical k values are ujq times a significant 
numerical factor. A simple dimensional analysis could 
have suggested /j, ~ k ^ luq. 

In Fig. [2] we show the continuum-case eigenfunctions 
Qui,u for uj = ±1, assuming fi = 50. As remarked earlier, 
Quj.u as a function of < u < 1 is a circular arc in the 
complex plane for fixed uj. We have normalized Quj,u so 
that the circle centers for different to lie on a vertical line 
in the complex plane and the function Qu]=o,u is a circle 
of unit radius. 

The example in Fig. [2] is close to the upper instability 
range, i.e., roughly where neutrinos would start convert- 
ing when the system evolves from high to low /i values 
in a SN. The angular modes remain fairly close to each 
other, corresponding to the observation that the multi- 
angle system evolves almost like the single-angle one. 

As time goes on and the unstable mode grows exponen- 
tially, the system is described by an eigenvector such as 
the one shown in Fig. [21 independently of the initial con- 
dition. The evolution consists of an exponential growth 
with rate k — Im(ri) away from the origin in Fig. [5] and 
at the same time a precession around the origin with fre- 
quency 7 = Re(ri). 



C. Normal- hierarchy solution 

We next repeat this exercise for NH. Here, in the single- 
angle case, the system is always stable, as already known 
from the flavor pendulum [8,]. However, in the multiangle 
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FIG. 2: Eigenfunction Qt^.n for a box spectrum at = 50. 
The solid arcs correspond to the physical range < it < 1. 
We also show as open circles an example for discrete angle 
eigenvectors with A'^^ = 20. 



case the system does have instabilities [T^] . In Fig. [T] we 
show the growth rate as a function of [i also for NH. The 
system is unstable only for a relatively small range 
and the instability parameter k is an order of magnitude 
smaller than for IH. 

Thus far we have ignored the ordinary matter effect. It 
can suppress self-induced flavor conversion [l7l - [l9| . and 
this "multiangle matter effect" is particularly effective in 
suppressing the NH instability. Therefore, we continue 
to focus on the IH case. 



IV. DISCRETE ANGLE MODES 

A. Eigenvalue equation 

We next turn to our main topic, the behavior of the 
system discretized in angle with Na > 1 angle bins, and in 
energy with Ne energy bins. As we always consider both 
neutrinos and antinautrinos, the number of frequency 
bins is — 2Ne- The spectrum is then implemented 
as 



s{uji - s{ub ■ 



i=l 6=1 



leading to 



i=l b=l 



(16) 



(17) 



One can then determine the eigenvalues f2 by solving 
Eq. ([TT]) which amounts to finding the roots of a polyno- 
mial in n of order N^Na- 

Alternatively, one can begin with the eigenvalue equa- 
tion of Eq. ([7]) in discrete form 

(Wfe + Uc\ - ^) Qk,c = E E + Si.b Qi.b ■ (18) 

i=l &=1 

This equation is of the form (M — il)Q = where Q is 
an N^Na dimensional vector of complex numbers and M 
a N^Na X NujNa matrix. What remains is to find the 
eigenvalues and eigenvectors Qn of M. Both methods 
provide the same eigenvalues and eigenvectors. 



B. Hair-comb spectrum 

We will concentrate on the discrete monochromatic 
spectrum (N^^ = 2) with Na angle modes representing 
the original box spectrum in the form 



1 ^" /h 
° fc=i 



1/2 



Na 



(19) 



For Na = I this is our previous single-angle case and 
for Na —7- oo we expect to recover the continuum limit. 
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We can solve the discrete version of the eigenvalue equa- 
tion in Mathematica without problem and we show the 
spectrum of growth rates in Fig. [3l 

We always find the "ordinary mode" which, for large 
Na, approaches the continuum solution (red dashed line 
in Fig. [3|). In addition, we find Na — 1 "extraordinary 
modes", which arise at larger /i. With increasing Na, 
the extraordinary modes shift their instability regions to 
larger ^ values and, in the limit Na — t- oo, disappear at 
infinity. Of course, for any finite Na, the extraordinary 
modes exist at a sufficiently large /i. 

In numerical studies, the neutrino flavor content is 
evolved along the radial SN direction, from large to small 
jj, values. A numerical integration must begin at a depth 
where the ordinary mode is stable, i.e., to the right of the 
red dashed curve in Fig. [3] Inevitably, the system first 
encounters the extraordinary instabilities, leading to spu- 
rious solutions. Therefore, for a chosen inner boundary 
radius ro, the number of angle modes Na must be large 
enough so that the extraordinary instabilities disappear 
to depths below it. If ro is chosen closer to the neutri- 
nosphere, implying a larger ^, the required Na is even 
larger. The most economical choice for the inner bound- 
ary radius vq is at the large-/z beginning of the ordinary 
instability region. 



C. Nature of the extraordinary modes 

The presence of more than one unstable mode is not 
surprising: solving Eq. ([TT|) for a discrete case leads to 
a polynomial in O of order 2Na and to equally many 
eigenvalues, some or all of which can have an imaginary 
part. However, empirically the extraordinary modes are 
quite distinct from the ordinary one. 

One difference can be gleaned from the general form 
of the eigenvector Q^j^u, Eq. ([8]), and especially the reso- 
nance denominator uj + uX — ^ — in. Varying u between 
and 1 over a continuous or discrete set of values may 
lead one to encounter a resonance of approximate width 
k/X. We illustrate this point in Fig. |4] where we show 
an example for the modulus of Quj,u as a function of u; 
the unshaded range corresponds to the physical range 
< u < 1. For the ordinary mode, the resonance lies 
in the unphysical range, implying that |(5i^,u| does not 
vary much as a function of u. In other words, all angle 
modes are close to each other and evolve similar to the 
single-angle case. Avoiding the resonance is possible if 
the precession frequency is negative 7 < 0, and this is 
indeed the case for the ordinary (or quasi-single-angle) 
mode. 

The extraordinary modes, on the other hand, always 
seem to have 7 > 0, i.e., they all precess in the opposite 
direction. The resonance falls in the physical u range and 
one or a few angle modes have IQi^.u! much larger than 
the others. In a plot of the circular arcs as in Fig. [21 
the ordinary mode for < m < 1 traces out a small part 
of the circle, but the extraordinary ones trace most of 



the circle. For our hair-comb spectra, the width of the 
resonance as a function of u has the approximate width 
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FIG. 3: Growth rates for unstable modes of the hair-comb 
spectrum in inverted hierarchy without matter (A = 0) and 
e — 1/2. The number of angle modes is Na ~ 1, 2, 5, 10 
and 20 (top to bottom). As a red dashed line we show the 
continuum case. The top panel corresponds to Fig. [T] 
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of the spacing of the discrete u modes; i.e., typicahy one 
or two angle modes are on resonance. The different ex- 
traordinary modes differ in the angle mode that is on 
resonance. By their very nature, these modes are not 
quasi-single-angle and it is unsurprising that, when they 
have grown beyond the linearized regime, they tend to 
cause kinematical decoherence among angle modes. 

We stress that, as far as the stability analysis is con- 
cerned, going beyond one energy bin (two frequencies) 
is not necessary. Since the equations of motion, Eq. Q , 
are linear in frequency to, an average frequency can repre- 
sent the whole spectrum. A multienergy treatment does 
provide, of course, more eigenvalues with respect to the 
monochromatic study, but no additional complex ones. 
The number of possible unstable solutions is related to 
the number of "spectral crossings" [2] . The additional so- 
lutions introduced by additional energy bins are purely 
real, and as such have no importance in the stability 
analysis. A multienergy treatment is necessary only in 
simulations if one is to resolve spectral features. 



D. Other cases of extraordinary modes 

Extraordinary modes are not limited to discrete angle 
distributions. For example, a u spectrum consisting of 
two boxes (instead of two delta spikes) has one ordinary 
and one extraordinary mode. Two boxes are equivalent 



to a single box with a gap, and the fi range where the ex- 
traordinary mode appears migrates to larger fj, values as 
the gap is chosen to be narrower. Likewise, if we consider 
two delta spikes with a separation Am, then the solution 
for Am = 1/2 corresponds to the upper panel of Fig. [3l 
Au = corresponds to the single-angle case, and for any 
other value one obtains two solutions, the extraordinary 
one migrating to higher /i values for decreasing Au. 

Generally it appears that "sharp features," notably 
steps, in the angle spectrum cause extraordinary modes. 
(Of course we mean the following: extraordinary modes 
with nonvanishing growth rate k > 0.) However, it 
happens only for the ascending steps (for increasing u), 
and not for the descending ones. For example, a de- 
scending staircase spectrum has only the ordinary (quasi- 
single-angle) mode, while an ascending one has as many 
extraordinary modes as steps minus 1. If the steps 
are somewhat smoothed, we still get the extraordinary 
modes. The location of the extraordinary mode on the fi 
axis depends on how narrow the spectral feature is, and 
the maximum k depends on the magnitude of the jump. 

A more mathematical classification of these observa- 
tions is not available at present. In a realistic SN situ- 
ation, the continuous angle spectrum is not a box, but 
typically a smoothly varying broad distribution. As such 
it should not have any extraordinary modes, or at least 
none with k values comparable to the ordinary mode. 
In realistic numerical SN simulations, the appearance of 
extraordinary modes is probably caused only by a dis- 
cretized angle spectrum. 
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FIG. 4: Modulus of the eigenvector Quj,u, according to Eq. ((8|, 
for fi — 50 and the hair-comb angle distribution with Na ~ 20. 
Only the range < m < 1 is physical. Top: Ordinary mode. 
Bottom: The extraordinary mode with the largest k. 



V. MULTIANGLE MATTER SUPPRESSION 

It is well known that a large density of matter, repre- 
sented by the parameter A, str ongly rn odifies the instabil- 
ity in a multiangle treatment [17h21| . This "multiangle 
matter effect" does not appear in the single-angle case 
and is one important motivation to go beyond single- 
angle studies in the first place. Typically, nonvanishing 
growth rates appear only for /x ~ A [la, [13 ■ In other 
words, for /Lt ^ A the instability is suppressed and collec- 
tive flavor conversion may not occur, notably during the 
accretion phase [13, [3 ■ 

To understand the impact of matter in the context of 
a discrete angle spectrum, we consider once more a box 
spectrum of angles. In Fig. [5] we show, as a red dashed 
line, the growth rate for different values of A. We see 
that indeed the low-/i instability is suppressed and that 
the instability region is shifted to fi ^ X. 

In addition, we show the solutions for a hair-comb rep- 
resentation with Na = 20. For the extraordinary modes, 
the presence of A has somewhat the opposite effect of en- 
hancing the growth rates at low /x values. In the bottom 
panel of Fig. [5] we see that the ordinary mode disappears 
in the forest of extraordinary ones. In other words, in 
the presence of matter one needs a yet larger value for 
Na to shift these modes away. 
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FIG. 5: Growth rates using the indicated values A for the 
matter effect. Red dashed linea: Continuous box spectrum 
for angles. Blue solid lines: Hair-comb representation with 
Na = 20. 



VI. DISCUSSION 

We have used the method of linearized stability analy- 
sis to shed light on the appearance of spurious solutions in 
numerical studies of nonlinear neutrino flavor evolution. 
Spurious solutions appear when the neutrino radiation 
field is represented by discrete modes and the number 
of angles Ng, is too small. The physical solution tends 
to be one where the angle modes remain nearly aligned 
(quasi single angle) and which we have called the ordi- 
nary mode. In addition, extraordinary modes appear at 
values of the effective neutrino density fi that depends on 
the spacing of the discrete modes. If they are sufficiently 
densely spaced, the extraordinary modes are at large fi 
(small SN radius) such that they do not spoil a numerical 
multiangle simulation. 

The linear stability approach allows one to determine, 
without much effort, the growth rate of the physical in- 
stability as a function of the SN radius. Any instability 
is important only if the growth rate is large enough to 
take the system into the nonlinear regime on the avail- 
able length scale. In a SN core, the effective neutrino 
and matter densities vary as power laws so that the rel- 
evant length scales correspond approximately to r. Typ- 
ical growth rates n are a few times the vacuum oscilla- 
tion frequency wq- For SN neutrino oscillations driven 
by the atmospheric mass difference, typically we have 
luq ^ 0.5 km^^ and so typical k values are few inverse 
kilometers. The available length scales are tens to hun- 
dreds of kilometers, so the physical instability has enough 
time to become nonlinear, in agreement with numerical 
studies. Of course, the number of e-foldings required for 
a mode to become nonlinear depends on its initial ampli- 
tude. 

In principle, then, for a concrete numerical example 
it is enough to find K(r) for the physical mode, based 
on a continuous angle distribution, and in this way find 
the onset radius of the instability. It would be enough 
to start the numerical integration at that radius. Since 
the starting point of flavor conversions is an exponen- 
tial runaway, nothing is gained by starting deeper, i.e., 
in the stable regime of the ordinary mode. The num- 
ber of angle modes then has to be chosen large enough 
that the growth rate of the extraordinary modes is much 
smaller than that of the ordinary one in in the onset re- 
gion. Starting the integration at a smaller radius requires 
enough angle modes to avoid the extraordinary modes 
becoming nonlinear before the system has reached the 
physical onset radius. 

If one determines the required Na by trial and error for 
a given numerical example, it has been observed that the 
solution becomes reproducible for Na above some criti- 
cal value, and nothing much changes by choosing Na yet 
larger. This behavior corresponds to the aforementioned 
requirement that the extraordinary modes must not be- 
come nonlinear before the physical onset radius. The 
solution then no longer changes because, at the onset ra- 
dius, the system will always select the physical eigenvec- 
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tor from whichever configuration the neutrino ensemble 
is in, by subjecting it to the exponential growth that will 
allow it to dominate the final outcome. 

The appearance of unphysical modes in the discretized 
problem suggests that one is using a bad representation 
of the physical system. It would be desirable to cast the 
original equations of motion in a form that allows for a 
numerical treatment while avoiding spurious solutions. 

Instead of using discrete angles one may consider 
spherical harmonics [l^]. We have attempted this ap- 
proach, truncating the expansion at some multipole or- 
der. Nevertheless, we obtain similar extraordinary modes 
and the required number of polynomials corresponds ap- 
proximately to the required number of discrete angles. 
It remains to be seen whether a suitable closure of the 
multipole equations of motion can be found that avoids 
unphysical solutions without going to very high multi- 
pole order. Conversely, if this is not possible, it would be 
important to understand whether extraordinary modes 
are unavoidable in any scheme that fails to resolve suffi- 



ciently fine details of the angle distribution. 

For now the question remains open if numerical brute 
force is the only way forward to understand SN neutrino 
flavor evolution in more general situations beyond the 
toy models that have been studied thus far, or if the 
equations can be set up in a way that avoids the need for 
excessive computer power. 
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